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We investigate the non-equilibrium relaxation dynamics of a one dimensional system of interacting 
spinless fermions near the XXZ integrable point. We observe two qualitatively different regimes: 
close to integrability and for low energies the relaxation proceeds in two steps (prethermalization 
scenario), while for large energies and/or away from integrability the dynamics develops in a single 
step. When the integrability breaking parameter is below a certain finite threshold and the energy 
of the system is sufficiently low the lifetime of the metastable states increases abruptly by several 
orders of magnitude, resembling the physics of glassy systems. This is reflected in a sudden jump 
in the relaxation timescales. We present results for finite but large systems and for large times 
compared to standard numerical methods. Our approach is based on the construction of equations 
of motion for one- and two-particle correlation functions using projection operator techniques. 
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Due to the peculiarities of quantum dynamics, it is 
possible to calculate expectation values of observables in 
the asymptotic state the system would reach after an infi¬ 
nite evolution without requiring the dynamics that lead, 
starting from a given non-equilibrium initial condition, to 
such state [1] . For non-integrable systems (even infinites¬ 
imally close to an integrable point) it has been shown that 
one-particle observables in such asymptotic state are in 
very good agreement with the predictions of statistical 
mechanics in the thermodynamic limit [2]. There is also 
evidence that the indicators of quantum chaos related to 
the statistical properties of the Hamiltonian’s eigenspec- 
trum [3] change abruptly when infinitesimally breaking 
integrability. However, thermalization of a closed system 
is a dynamical process. In particular, the dynamics to¬ 
ward the final thermal state in nearly integrable systems 
seems to be more rich and subtle than the dichotomic 
situation present in the static results. 

In fact, it is well known that systems close to an inte¬ 
grable point may exhibit non-thermal stationary states, 
both in one dimension (ID) and in higher dimensions [4- 
12]. Such metastable states emerge as a result of a fast 
lost of memory of the initial conditions due to dephas¬ 
ing [11, 13]. The emergence of such non-thermal sta¬ 
tionary states has been observed experimentally in cold 
atomic Bose gases [14, 15]. Quite interestingly, dephas¬ 
ing alone may lead to complete thermalization of some 
observables, such as the kinetic energy of the system, 
a phenomenon dubbed prethermalization in Ref. [16]. 
However, little is known about the subsequent evolution 
of the system after getting caught in such metastable 
states, in part due to the formidable technical challenge 
that poses the out-of-equilibrium dynamics of many-body 
quantum systems. For example, simple questions such as 
how much do these stationary states live or which are the 
relevant timescales involved in the thermalization process 
has no definitive answer yet. 


In this work we analyze the dynamics of a system of 
ID spinless fermions near the XXZ integrable point while 
varying the distance to integrability and the energy of 
the system. The main finding of our work is that be¬ 
low a certain finite threshold distance to integrability 
and for sufhciently low energies these metastable states 
are extremely long-lived. Their lifetime increases by sev¬ 
eral orders of magnitude while crossing such threshold. 
Such metastable states may completely hinder the obser¬ 
vation of the final thermal equilibrium state in experi¬ 
ments or numerical simulations leading to an apparent 
lack of thermalization. We find the situation rather sim¬ 
ilar to that of glasses, systems that exhibit the typical 
two-step relaxation as a consequence of getting caught 
in extremely long-lived metastable states, whose lifetime 
can be of geological scale for sufficiently low tempera¬ 
tures or high densities. Glass-like behavior, including 
ageing and slow relaxation, has also been found in open 
quantum systems [17, 18]. Our approach, based on the 
projection operator formalism [19-24], enables to inves¬ 
tigate not only the dynamics of the off-equilibrium mo¬ 
mentum distribution, but also of two-times correlation 
functions for large times and system sizes. This allows to 
envisage the rich relaxation dynamics of nearly-integrable 
systems, which turns out to be characterized by several 
relevant timescales that we quantitatively study in the 
specific example at hand. 

We consider a ID model of spinless fermions with 
nearest neighbor interactions, and nearest and next-to- 
nearest neighbor hopping, H{Ji, J 2 , A) = Ho{Ji, J 2 ) + 
Hi{A), 

L-1 L-1 

Do(Ji, J 2 ) = -Ji {c]cj+i+h.c.)-J 2 (c^Cj+2+h.C.), 

j=0 j-0 
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^i(A) = A^ 

(^) 

i=o 

where L is the number of sites in the chain, the c op¬ 
erators obey canonical anticommutation relations and 
rij = CjCj. We assume periodic boundary conditions, 

^L+m = ^2 = 0 the model is integrable through 

Bethe ansatz, while for any J 2 ^ 0 the model is non- 
integrable. 

We are interested in the dynamics of the system 
starting from a non-equilibrium initial condition p( 0 ), 
[H,p{0)] ^ 0. We shall restrict to the weakly inter¬ 
acting case a = A/Ji 1 and consider only ho¬ 
mogeneous initial states. In such case the operators 
defining the momentum space density distribution = 
r j ^ = 0 ,..., L — 1 , emerge as 

the natural slow variables of the system. Using the pro¬ 
jection operator technique it is possible to derive an ex¬ 
act coupled system of equations of motion for the av¬ 
erage n{k,t) = {nk{t)), where (...) = Tr[p(0)...], and 
the fluctuations F{k,k',t) = {nk{t)nk') — {fik(t)){nk>) 
of the slow variables [19-24]. We work in the Heisen¬ 
berg picture, 0{t) = (^ = !)• The integro- 

differential equations are manageable if (i) we restrict 
to uncorrelated initial conditions of the form p( 0 ) = 
with Z = Tr[e-^'=^('=’°)”''], which in¬ 
clude free fermion initial states, and (ii) we keep only 
the leading order in a. In particular, the equation for 
the averages n{k, t) is formally identical to that derived 
using heuristic arguments in Ref. [10]. For the explicit 
expressions and a thorough derivation of such evolution 
equations see [24, 25]. Although the approach involves 
perturbative steps, the results go beyond conventional 
lowest order perturbation theory because the perturba¬ 
tion expansion is performed inside the integro-differential 
equations and, therefore, the coupling is involved in a 
highly non-linear way in the final expressions. The equa¬ 
tions of motion can be efhciently solved numerically al¬ 
lowing to study large systems (L ~ 10 ^) and times far 
beyond the reach of standard numerical techniques such 
as time-dependent density-matrix renormalization group 
(t-DMRG), which allows to access the long time dynam¬ 
ics of the system in the thermodynamic limit. We fix 
a = 0.2 for the rest of the paper. 

We begin by considering the effect of varying the 
strength of the integrability breaking parameter J 2 . 
As initial state we consider the ground state j'I'o) 
of Hq{Ji,Q) with momentum distribution n{k,0) = 

0{eo{k) — eoikp)), where eo(fc) = —2Ji cos( 2 A: 7 r/L). We 
will work at half-filling, kp = L/A. We shall first fo¬ 
cus on the relaxation of the momentum modes close to 
the Fermi surface, whose relaxation timescales are the 
longest among all momentum modes. In particular, we 
will study the decay of the quasiparticle residue, defined 
as^(t) = n(kF,t)—n{kp-\-l,t). In higher dimensions this 
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FIG. 1. Decay of the quasiparticle residue Z{t) varying J 2 
for a system with L = 1024 away from any revival. Upper 
panel: Plot in semilogarithmic scale. Dashed (red) lines are 
exponential fits. Lower panel: Semi logarithmic plot showing 
the details of the decay for J 2 < 0.5Ji, the curve for J 2 = 
O. 7 J 1 is also presented for comparison. Dashed (red) lines 
are stretched exponential fits. Inset: Log-log plot. Dashed 
(black) lines are power laws with the exponent given by the 
LM predictions Z{t) ~ t with 7 = K~^ — 2). 

The Luttinger parameter is obtained from bosonization, K = 

where wf is the Fermi velocity. 


quantity exhibits typical “prethermalization plateaus”, 
and has become a standard tool to detect metastable 
states in fermionic systems [4, 5, 11]. In ID this quan¬ 
tity does not exhibit plateaus [26]. We shall show that, 
nevertheless, metastable states are present and that they 
profoundly affect the relaxation of the momentum modes 
close to the Fermi points. In the upper panel of Figure 1 
we show the decay of the quasiparticle residue for dif¬ 
ferent values of J 2 for a system with L = 1024 and for 
times far away from any recurrence effect. In such sit¬ 
uation finite size corrections are negligible and we are 
thus accessing the thermodynamic limit dynamics. It is 
clearly visible that for J 2 < 0.5Ji the decay of Z{t) is 
extremely slow (subexponential) while for J 2 > 0.5Ji it 
exhibits exponential behavior. In the lower panel we an¬ 
alyze in more detail the slow evolution. In particular, 
we find that after a fast Gaussian-like initial evolution 
taking place for tJi < 1 , there is a faster than power law 
but slower than exponential decay. In this regime it is 















3 


possible to fit a stretched exponential Z{t) = 
where both r and j5 depend on A and J 2 . In the fits 
showed in the lower panel of Figure 1, r Ji ~ 10^ and 
(3 ~ 0.27. On the other hand, for the exponential de¬ 
cays taking place for J 2 > 0.5Ji, tJi ^ 10^. Below we 
shall see that such abrupt decrease of the relaxation time 
scales is related to an abrupt increase in the lifetime of the 
metastable states. We finally note that for J 2 < 0.5Ji, 
after the initial Gaussian evolution, the initial trend of 
the curves is very well described by the power law decay 
predicted by the non-equilibrium dynamics of the Lut- 
tinger model [27, 28] . In the t-DMRG study in Ref. [29] , 
where this remarkable fact was first noticed, such initial 
trend was the only accessible portion of the dynamics. 
The fact that the evolution equations capture such fea¬ 
ture of the dynamics clearly indicates that the results are 
not perturbative in the usual sense. 

Two-time correlations turn out to be an adequate tool 
to detect and analyze the prethermalized states. In par¬ 
ticular, we introduce the connected correlation function 
of the current operator J = Yi ~ 

Cj{t) = {Jit)J)-{J{t)){J), (3) 

which can be obtained from the fluc¬ 
tuations of the slow variables Cj{t) = 

sin(2fc7r/L) sin(2fc'7r/L)F(fc,/c', f). It is im¬ 
portant to note that two-times connected correlation 
functions of local observables are a standard tool to 
diagnose the ergodic status of equilibrium dynamics, 
both in quantum [30-32] and classical systems [33]. 
Loosely speaking, this type of correlators is expected 
to decay rapidly to zero for ergodic systems, since the 
initial state and the state at time t are expected to 
be completely decorrelated after some characteristic 
timescale, i.e., the system is expected to loose memory 
of the initial condition. For non-ergodic systems the 
correlator is expected to saturate to a non-zero constant. 
An intermediate behavior, with a plateau emerging in 
between a first fast evolution and the final decay to 
zero, arises in systems that get caught into long-lived 
metastable states, the most prominent example being 
glassy systems, such as spin glasses and supercooled 
liquids [33, 34]. On general grounds we expect that 
this kind of correlators should perform as a similar 
diagnosis tool in the non-equilibrium situation under 
consideration. In the left panel of Fig. 2 we show the 
decay of Cj{t) for the ground state initial condition 
varying the value of the integrability breaking parameter 
J 2 in a system with L = 800 and for times away from any 
recurrence effects, which, again, amounts to investigate 
the dynamics in the thermodynamic limit. For large 
values J 2 Ji the decay develops in a single step. 
Decreasing J 2 a plateau arises in between the initial 
Gaussian evolution and the final decay. The length of 
the plateau becomes larger as we further decrease J 2 yet 
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FIG. 2. Decay of the correlation function Cj{t). Left panel: 
starting from the T = 0 initial condition varying J 2 , for which 
the energy density is e = 0 independently of the value of J 2 . 
Inset: Relaxation of the kinetic energy for the same param¬ 
eters as the main figure. Right panel: for fixed J 2 = 0.2Ji 
varying the temperature of the initial state. Note the loga¬ 
rithmic scale on the time axis. 



it does not increase smoothly but rather seems to have 
an abrupt jump exactly for J 2 = 0.5Ji. This clearly 
indicates that the system is caught in non-thermal 
metastable states whose lifetime grows abruptly around 
J 2 ~ 0.5 Ji. 

Another factor that deeply influences the lifetime of 
such prethermalized states turns out to be the energy of 
the system. In particular, we investigated the behavior 
of Cj{t) starting from finite temperature initial states 
n(fc, 0) = (1 -I-exp[(e(fc) — e{kF))/T])~^, where T is the 
temperature (ks = 1). The energy density of the system 
e = rTr[p(0)F]-eo, with cq = ^ (Toli7( Ji, 0, AjjTo), is 
a smooth, monotonous function of T. In the right panel 
of Fig. 2 we show the relaxation of Cj{t) for J 2 = 0.2 
varying the temperature of the initial condition. The ef¬ 
fect of increasing the temperature (energy) of the initial 
state is to gradually decrease the lifetime of the prether¬ 
malized states. In particular, for sufficiently high ener¬ 
gies, they are completely suppressed. We observe that 
the metastable states emerge in the same timescale in 
which the kinetic energy of the system ekin{t) = {Ho{t)) 
saturates to its final value (see inset in Fig. 2). Such 
prethermalization timescale turns out to be independent 
of the value of J 2 and the energy of the system: tptJi ^ 5 
in all cases. 

A simple qualitative picture can be formulated. The 
initial fast evolution taking place for tJi < 1 is caused 
by the dephasing of some quasifree modes of the sys¬ 
tem [35]. For fermionic systems these can be identified 
with the bosonic modes associated with the bosonization 
of the excitations close to the Fermi points surface [11]. 
This initial regime is not sensitive to the details of the 
interaction and is analogous to the initial ballistic ex¬ 
pansion in the relaxation of classical glasses and fluids 
in general. The subsequent relaxation of the system is 
provided by inelastic collisions. If the inelastic relax¬ 
ation channels are scarce the system becomes trapped in 
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FIG. 4. Characteristic decay timescales of the correlator 
Cj{t) for three different representative values of J 2 as a func¬ 
tion of the energy of the initial state. Inset: Decay timescales 
of the dynamical distance di 2 as a function of the energy den¬ 
sity of the pair of initial conditions for three representative 
values of J 2 , 0, 0.2Ji and O. 8 J 1 (with the same symbols as 
the main figure but with dashed lines). 


FIG. 3. Top left panel: Decay of the dynamical distance di 2 (t) 
(note the logarithmic scale in the time axis) for J 2 = 0.2Ji 
varying the energy of the pair of initial conditions. Top right: 
Overall decay of the correlation function Cj{t) for J 2 = 0.2Ji 
varying the temperature T of the initial state. Bottom panel: 
Detail of the behavior for Cj{t)/Cj{G) ~ 1 making visi¬ 
ble the presence of two-step relaxation and long-lived quasi¬ 
stationary states (note the logarithmic scale on the time axis). 


metastable states whose lifetime is a measure of the rate 
of occurrence of such inelastic scattering events. For the 
model H{Ji, J 2 , A) it can be shown that as soon as J 2 
becomes larger than 0.5 Ji the manifold of kinetically al¬ 
lowed collisions (those that conserve momentum and ki¬ 
netic energy) is dramatically enlarged [36]. Such type of 
collisions are included in (but do not exhaust) the equa¬ 
tions of motion that we consider [25]. The number of 
relaxation channels is thus drastically enlarged beyond a 
finite threshold away from integrability. Moreover, since 
at higher energies there are more possible inelastic colli¬ 
sions, the lifetime of the metastable states is suppressed 
as we increase the energy of the system. 

In order to make quantitative statements about the 
lifetime of the prethermalized states we find convenient 
to concentrate on smaller systems, with L = 256, which 
shall allow us to access longer times. In this case it 
must be noted that, specially for long times, deviations 
from the thermodynamic limit dynamics may be appre¬ 
ciable [25]. To characterize the overall relaxation of the 
momentum distribution we study the dynamics of two 
replicas of the system. In particular, we prepare two ini¬ 
tial conditions ni(fc,0) and n 2 (fc, 0 ) with approximately 
the same energy and particle density, and define the dy¬ 
namical distance between the time evolved momentum 


distributions 

(4) 

k 

For a system whose only conserved quantities are total 
energy and particle number this distance should decay 
to zero at long times. For systems with additional con¬ 
servation laws (like the integrable model J 2 = 0) it may 
saturate to a non-zero constant. In any case, the presence 
of short time plateaus in the time evolution of d\ 2 {t) rep¬ 
resents a clear sign of the formation of metastable states. 
In other words, if the two systems, prepared in differ¬ 
ent initial conditions, get caught in different metastable 
states that are at a distance d\ 2 ^piat of each other, then 
a plateau in (ii 2 (t) at the value piat would be present. 
If plat 0 we can be sure that these are non-thermal 
metastable states, at least for the non-integrable model. 
We prepare the initial conditions slightly perturbing free 
fermion thermal states. In the top left panel of Fig. 3 
we show the decay of ^ 12(0 for L = 256 and J 2 = 0.2Ji, 
but the results are similar for any J 2 < 0.5Ji [37]. We 
find that if the energy of the pair of initial conditions 
is low enough di 2 {t) relaxes in two steps. The length of 
the plateau increases for lower energies. In contrast, for 
J 2 > 0.5 Ji we find single-step relaxation in all the energy 
range (not shown), confirming the picture that emerged 
from the analysis of Cj{t). In the top right panel of Fig. 3 
we show the overall relaxation of the correlation function 
Cj{f) also for J 2 — 0.2Ji. For high energies we find that 
it decays to zero on the accessible timescales. For low en¬ 
ergies we observe a very pronounced slowing down of the 
relaxation. In the bottom panel we show that the slowing 
down is caused by the presence of long-lived metastable 
states whose large lifetime can be appreciated. 































5 


In Fig. 4 we show the decay timescales of Cj{t) for 
different values of J 2 . The relaxation timescale r was de¬ 
fined as C'j(t)/Cj(0) = 0.6, in order to extract the max¬ 
imum possible number of data points from the results 
at disposal while still considering a faithful indicator of 
the relaxation timescale. We see that for J 2 < 0.5Ji the 
relaxation timescales show an abrupt increase for suf¬ 
ficiently low energies. For J 2 > 0.5Ji the relaxation 
timescale is almost unchanged as we vary the energy of 
the initial condition. 

A special remark is in order with respect to the effec¬ 
tive decay timescales of d\ 2 it)- Being quite independent 
of the specific form of the initial conditions, this is one of 
the characteristic timescales in the thermalization pro¬ 
cess [16]. It is the timescale beyond which the system 
has lost all memory of the initial conditions at the level 
of single particle observables. Nevertheless this does not 
mean that the system is in equilibrium. In fact, we find 
that t[Cj\ is, at least, one order of magnitude larger than 
T[(ii 2 ] (defined in the same way as for Cj) in all cases in 
which we have data to compare. This is illustrated in the 
inset of Fig. 4. The fact that one-particle observables 
attain thermal behavior much before than two-particle 
correlations indicates that thermalization is a hierarchi¬ 
cal process. The information about the initial conditions 
encoded in a correlation function increases with its or¬ 
der. Our results suggest that, accordingly, the relaxation 
timescales also increase monotonically with the order of 
the correlation function. However, we may also expect 
that for some finite (but possibly very large) order corre¬ 
lations do not thermalize at all, reflecting the unitarity of 
quantum dynamics. Finally, we note that it is for times 
larger than T[(ii 2 ] that a description based on kinetic 
(memoryless) equations, such as the quantum Boltzmann 
equation [38], is justified [39]. 

This work was partially supported by CONICET (PIP 
0662), ANPCyT (PICT 2010-1907) and UNLP (PID 
X497), Argentina. 
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Supplementary Material to: Glass-like Behavior in a System of One 
Dimensional Fermions after a Quantum Quench 

In this supplement we will show the derivation of the evolution equations used to extract the results discussed in 
the main text. We will make a brief outline of the derivation of the evolution equation for the momentum distribution 
since it has been already been presented in full detail in Ref. [24], and pay most of the attention to the evolution 
equation for the fluctuations. We also include a discussion on the behavior of the relaxation timescales with system 
size. 


EVOLUTION EQUATION FOR THE MOMENTUM DISTRIBUTION 


We shall first make some elemental definitions to set up the situation. We consider a system of interacting spinless 
fermions with Hamiltonian H = Hq + aHi with 

H = Ho + aHi ='^e{k)n{k) + a ^ V^^'l^^c\ki)c\k 2 )c(ko)c{k 4 ), (1) 


where c^(fc) and c{k) are fermionic creation and annihilation operators satisfying canonical anticommutation relations, 
is the momentum-space matrix element of the interaction, e(fc) is the dispersion relation, n{k) = {k)c{k) is the 

number operator and a is the strength of the interaction. Our results can be easily extended to the bosonic case. The 
hermiticity of the Hamiltonian and the symmetry in the sum indices impose Vk^k^ = — Vk^k^ = — Vk^k^ = Vk^k^, 
where V denotes the complex conjugate. 

Furthermore, we will be interested in the special case of a translationally invariant Hamiltonian in which the 
particles interact via a pair potential v(x — y). In such case 

= ^'^fci+fc2.fc3+fc4 (^(fei - ki) - vik2 - ki) - v{ki - fca) -b v{k2 - fca)), (2) 

where v{k) is the Fourier transform of the potential and we have written the antisymmetrized version in order to 
respect the symmetry conditions of the potential. We are interested in the evolution of the system starting from an 
arbitrary initial condition given by a density matrix p(0) which we leave unspecified for the moment. 

In the specific ID model treated in the main text e{k) = —2 Ji cos(2/c7r/L) — 2 J 2 cos(4fc7r/L) and v{k — k') = 

To obtain an evolution equation for the momentum distribution we start from the Liouville equation in the inter¬ 
action representation {h=l): 

dtpit) = -ia[Hi{t),p{t)] = aL{t)p{t), (3) 

where 0{t) = is the interaction representation of the operator O and we have introduced the Liouville 

superoperator L{t)0 = —i[iLi(t), O]. Our task is to find approximate solutions to the microscopic dynamics described 
by the Liouville equation. The POT defines a program for achieving this. We need to first identify the “slow” or 
“macroscopic” variables in our system and then project the dynamics into the subspace of these slow variables. As 
noticed in the main text, in a weakly interacting homogeneous system the occupation number operators emerge as 
natural slow variables since [H,n(k)] = 0{a). To perform the projection we first introduce the “relevant” density 
matrix 


a{t) 




'^X{k,t)n{k) , 

k 


( 4 ) 


where the time-dependent partition function is given by Z{t) = Tr [exp A(A;, t)n(fc))]. Note that a{t) = cr{t). 

The Lagrange multipliers X{k,t) enforce the relation: 


(n(fe))t = Tr[n(fc)(T(<)] = Tr[n(fc)p(t)]. 


( 5 ) 


The projection of the dynamics consists in finding an equation of motion for a(t). To this end we introduce a 
projection super-operator P{t) that projects the relevant density matrix P{t)p{t) = d{t): 


P{t)p = ( a{t) - ^ 


Sa{t) 

5{n{k))t 


r(fe))* TrM + ^ 


S(7{t) 

6{n{k)) 


-Tr [n{k)p] 


( 6 ) 
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where /i is an arbitrary density matrix. The projection operator (6) is specially designed to satisfy the following 
properties [19-21]: 


P{t)p[t) = a{t), 

P{t)dtp{t) = dta{t), 

Tr [n[k)P{t)p\ = Tr [n{k)p\ , 

Pit)P{t')p = P{t)p, (7) 

P{t)L{t)P{s)p = 0. (8) 

The fourth identity, setting t = t', expresses the idempotent character of the projector, while the last identity depends 
on the explicit form of the Hamiltonian 77, in particular, on momentum conservation. It is also useful to define the 
complementary projector Q{t) = 1 — P(t). 

Following the usual steps [19-22], introducing projectors in the Liouville equation, we obtain an equation for the 
dynamics of the slow degrees of freedom 


dtP{t)p{t) = aP{t)L{t)p{t), (9) 

and other for the fast, microscopic degrees of freedom 

dtQ(t)p{t) = aQ{t)L{t)p{t). (10) 

Inserting the identity I = P{t) + Q{t) in both equations we obtain the system: 

dtP(t)p{t) = aP{t)L{t)P{t)p(t) + aP{t)L{t)Q{t)pit), (11) 

dtQ{t)p{t) = aQ{t)L{t)P{t)p{t) + aQ{t)L{t)Q{t)p{t). (12) 


The equation for the relevant density matrix a{t) can be obtained solving the equation for the irrelevant part Q{t)p{t) 
in the second line of the system and inserting the solution in the first line. The second line is a linear first order 
homogeneous differential equation in the operator Q{t)p{t) (the inhomogeneity is aQ{t)L{t)P{t)p{t)) that can be 
(formally) solved in the same way as a real valued function differential equation. The solution is: 

Q{t)p{t)=af dsG{t,s)Q(s)L{s)P{s)p{s) + G{0,t)Q{0)p{0), (13) 

Jo 


where G{t, s) is an ordered exponential G(s, t) = T_,. exp 


—a f* ds' Q{s')L{s') 


i.e., the solution of the equation 


dtG{s,t) =-aG{s,t)Q{t)L{t), (14) 

G{s,s) = I. (15) 


Inserting (13) in the first line of (11) we obtain the desired equation: 


dtd(t) = aP{t)L{t)d{t) + f ds P{t)L{t)G{t, s)Q{s)L{s)d{s) + aP{t)L{t)G{t,0)Q{0)p{0), (16) 

Jo 

The first term in Eq. (16) is a mean field-like term that vanish due to momentum conservation, the second one can be 
expressed entirely in terms of the past history of the momentum distribution (n(fe))t, and the third one is a microscopic 
noise that can not be expressed in terms of the slow variables. The last term in Eq. (16) (the microscopic noise) 
disappears if we chose an initial condition of the same form of the relevant density matrix, i.e., if p(0) = cr(0). We shall 
then chose Gaussian (uncorrelated) initial density matrices, such as the ground state of Hq or a finite temperature 
state. We are thus considering an interaction quench. 

Eq. (16) is equivalent to the Liouville dynamics and, in general, as difficult to solve as the original problem. It sets, 
however, a good starting point for approximations. To render Eq. (16) tractable we perform a perturbative expansion 
in the interaction strength using that G{t, s) = I + 0{a). Taking the trace {n{k))t = Tr[n(A;)tT(t)] we finally obtain 

dt{n(k))t = f dsTr [n{k)L{t)L{s)a{s)] + 0{a^). 

Jo 


(17) 
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A great simplification arises since, given the Gaussian structure of cr(t), we can use the Wick pairing rule to evaluate 
the trace in (17). After a straightforward (but potentially tedious) calculation we obtain the explicit equation of 
motion 


I- I- ^0 


sin 




Ae 


k,ko 


/C2,fc3,fc4 ^ ‘~^'^k3,k4 

X {f{k, s)f{k 2 ,s)f{k 3 , s)/(fc 4 , s) - /(fca, s)f{k 4 , s)f{k, s)f{k 2 ,s)) + 0{a^), 


(18) 


where = e(fe) + £(^ 2 ) — £(^ 3 ) — £(^ 4 ) and, in order to ease the notation, we have defined f{k,t) = {n{k))t 

and f{k,t) = 1 — {n(k))t. This equation, in slightly different versions, has appeared many times in the literature. In 
Ref. [20] it was derived using the same tools that we present here [41] but it was used only as an intermediate step to 
derive the Boltzmann equation whereas in Refs. [10, 38] it was heuristically derived and used to study the dynamics of 
infinite dimensional models and to derive a quantum version of the Boltzmann equation, respectively. Eq. (18) is valid 
for systems in the continuum limit and also for lattice systems which only conserve quasi-momentum. A discussion 
on the accuracy of Eq. (18) can be found in Ref. [24]. 

It is worth noticing that Eq. (18) includes kinetic collisions = 0) as well as non-kinetic processes. This is 

related with the fact that Eq. (18) describes the dynamics of the system in all timescales. For short timescales, where 
there has not yet elapsed enough time for particles to collide, non-kinetic processes dominate the dynamics, whereas 
for long timescales kinetic collisions are dominant. 

With respect to implementation details, Eq. (18) can be solved using standard techniques for systems of Volterra 
integral equations [42] . A straightforward algorithm for the solution using, for instance, the trapezoidal rule to perform 
the time integral, implies a calculation time that scales as x where N is the number of times steps and D 
the space dimension. We have found an algorithm whose execution time scales as x N allowing us to reach large 
sizes and times. We finally note that the evolution equations are very suitable for parallel computing. 


EVOLUTION EQUATION FOR HIGHER ORDER CORRELATIONS 

In this section we undertake the calculation of two-times correlation functions of the slow variables, often referred 
as the fluctuations of the variable. We will briefly review the projection operator formalism of [19]. To perform the 
calculation we need to switch to the Heisenberg representation, where an operator with no explicit tome-dependence 
evolves according to the law: 


dtO = i[H,0] = iLO. (19) 

Note the difference with the Liouvillian defined in (3). The trace operation defines a dual projection operator over 
the observables of the Hilbert space: 


TT[OPit)^l]=Tl[^iP{t)0], 


( 20 ) 


where O is an observable, /i a density matrix and P(t) is the observable space projection operator. Its explicit form 
can be obtained right from the last expression and reads in our case. 


P{t)0 = Tr[<j{t)0] + ^(n(fc) - (n(fe))t)Tr 

k 


6a{t) 

_S{n{k))t 


O 


It can be readily shown that this dual projector satisfies the convenient properties: 


P(t)P(t') = P(t'), 

P(t) = P(t)P(t)(l-P(t)). 


( 21 ) 


( 22 ) 

(23) 


It is also useful to define the complementary projector Q{t) = 1 — P(t). 

In the Heisenberg representation the strategy is to separate the slow and fast components of the evolution operator 
using the projector P(t) and its complement Q(t): 


giLt 


= P^^P{t) + P^*Q{t). 


(24) 
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Using the last identity in Eq. (22) we obtain the equation for the irrelevant part 

dte^^^Qit) - e^^*Q{t)iLQ{t) = \iL - P(t)j Q(t), 

that can be solved formally using the propagator G{s,t) satisfying 

dtG{s,t) = -iLQ{t)G{s,t), 

G{s,s) = I, 


i.e., the chronologically ordered exponential 


G{s,t) = exp 


— / ds' iLQ{s') 


The solution is 


= e*-^"Q(s)G(t,s) + / dwe*'^“P(w) iL-P{u) Q{u)G{t,u). 


Inserting the formal solution into the relevant part we obtain 


=e*^‘P(t)+ / due^^^P{u){iL-P{u)){l-P{u))G{t,u) 

J S 

+e*^^(l-P(s))G(t,s), 

where s is an arbitrary time 0 < s < t. Defining a new ordered exponential 


(25) 

(26) 

(27) 

(28) 

(29) 


(30) 

(31) 


G(s, t) = r<_ exp 


ds' iLQ{s') 


(32) 


we can write 


=e*^‘P(t)+ / due^^''Piu){iL-P{u)){l-P{u))G{u,t) 
J S 

+e*i«(l_ P(s))G(s,t). 


(33) 

(34) 


Using the explicit form of the projector superoperator Eq. 21, it is possible to obtain an operator Langevin-like 
equation for the slow variables [19], 


n{k,t) = e*^‘?i(fc)= Vkit) + ^ Glk.k'{t) 6n{k,t) + 

k' 

+du ^Kk{t, m) + ^ (t, u) dn{k', u) 

+??/c(f,s). 


We have defined 


n{k) = iLn{k), 

6n{k,t) = n{k,t) — {n{k))t, 


(35) 

(36) 

(37) 


(38) 

(39) 


where n{k,t) = e^^*n{k) . The organized drift Vk{t) = Ti {n{k)[H, a(t)]} and the collective frequencies Glk.k'{t) = 
Tr s{n{k))t = s\nlk))t ’ vanish identically due to momentum conservation in Hi. The after effect functions 

Kk{t,u) can be expressed in terms of the {n{k))t’s: 


Kk{ii u) = Tr [a{u)iLQ{u)G{u, t)h{k)] = —a^Tr {[i7i(M), [Hi{t),n{k)]]a{u)} + 0{a^), 


(40) 
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and are related with the dynamics of the momentum distribution, see Eq. (17). In the last equality of Eq. (40) we 
have used G(u,t) = I + 0{a). The memory functions 4)fc,fc/(M,t) read 


= Tr 


Sa(u) 


6{n{k')) 
-^(n(fc"))«Tr 


iLQ{u)G{u, t)n{k) 


_S{n{k'))uS{n{k'')), 


-G{u, t)n{k) 


(41) 

(42) 


This functions are related to the after effect functions via a functional derivative [19] 


^k,k'{t,u) 


6 JgdsKk{t,s) 
8{n{k'))u 


(43) 


which constitutes the key relation between the dynamics of the momentum distribution and the fluctuations. Lastly, 
we have defined the microscopic noise 


77 fc(t, s) = e*-^^(l - P(s))G(s, t)n(k). 


(44) 


If we take the trace with respect to /^(O) in the Eq. (35) and we keep only with the lowest order in a we will recover 
the kinetic equation for the momentum distribution Eq. (18). But our intention is to rederive that equation but 
to make approximations on the dynamics of the operators themselves in order to calculate higher order correlation 
functions. 

Setting s = 0 in the Langevin Eq. (35) and subtracting the mean value we obtain an equation for the fluctuations: 


dtSn{k, t) 


where Sn{k,t) = e^^^n{k) — {n{k))t and 



'^^k,k'{s,t) 6n{k',s) +r]k{t), 

k' 


(45) 


rik{t) =?7fe(t,0)-Tr[p(0)?7fe(t,0)] =?7fc(t,0). 


(46) 


From this definition of the noise is clear that {r]k{t)} = Tr [p{0)r]kit)] = 0. 

Starting from Eq. (45), using the Wick rule and taking the functional derivative in Eq. (40) we find an explicit 
evolution equation for the time-correlation function of the slow variables Fk,k'{t) = Tr [p{0)Sn{k,t)6n{k',0)]: 


Fk,k'{t) = Fk,k'{0) — 16q!^ 



Fk,k'{s)Ak{t,s) + ^Eq,fe/(s) [Bk,q{t,s) - 2Ck,qit,s)] 


(47) 


where the matrices can be written as (using the notation defined earlier) 


Bk.qit,s)= KlJ 

k3,ki 

Ck.q{t^s)= Y, 

k2 ,k4 


Sin 




sin 


\t-s)^e^qt 


^ kM 


(/(fc, s)/(fe 3 , s)/(fe4, s) + f{k, s)/(fe 3 , s)f{ki, s)) , 


(/(fc, s)f{k 2 ,s)f{k 4 , s) + f{k, s)/(fe 2 , s)f{k 4 , s)) 


and Ak{t,s) = J2k' Bk',k{t,s). To arrive to Eq. (47) we have to take into account that ^k,k'{t,s) = 0{a'^) and that 
G{s,t) = I + 0{a). Notice that the matrices satisfy the convenient property Ak{t,t) = Bk.q{t,t) = Ck.q{t,t) = 0 
reflecting causality, i.e., the value of the fluctuations at time t only depends on the history of the fluctuations for 
times strictly before t. A similar statement can be done for the equation for the momentum distribution Eq. (18), 
the kernel of the integral equation vanishes for s = t. This property of the evolution equations brings a big technical 
simplification since it is not necessary to solve autoconsistent equations at each time step in the numerical integration. 

To the best of our knowledge, the equation (47) was first presented in Ref. [24] and used to investigate the dynamics 
of a concrete system in the present publication. In order to solve Eq. (47) we need first to know the dynamics of 
the momentum distribution, i.e., calculate the solution to Eq. (18), in order to determine the coefficients Afc(t, s), 
Bk,q{t, s) and Ck.q{t, s). With this input, Eq. (47) is as amenable to numerical solution as Eq. (18). We finally remark 
that since the projection operator in the Heisenberg representation works directly on the evolution operator itself it 
would be possible to obtain similar evolution equations for other observables. 
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FIG. 1. Relaxation timescale of the correlator Cj{t) as a function of the size of the system. Left panel: For J 2 = 0.2Ji and 
T = O. 3 J 1 , deep in the parameter region with long lived prethermalized states. Right panel: For J 2 = O. 8 J 1 and T = 0.5Ji, 
deep in the parameter region exhibiting one-step relaxation. 




RELAXATION TIMESCALES AND SYSTEM SIZE 

Having considered finite size results it is important to know about the dependence of the results on the system 
size, L. In Fig. 1 we show the dependence of the relaxation timescale of the correlation function Cj{t) with system 
size for a system deep in the “glassy” phase and for a system with normal relaxation. We observe that the timescales 
of the normal system increase almost linearly with a slope ^ 2 until, around L ^ 500 it begins to saturate to 
the thermodynamic limit value. For the glassy system the increase is again almost linear with system size with a 
(considerably larger) slope ~ 20 until, around L ^ 700 it begins to saturate to the thermodynamic limit value. This 
is another sign pointing to the fact that the relaxation mechanism in the glassy phase is qualitatively different from 
that of the normal phase. 




